Heat and Trees

Right click to download this notebook from GitHub.


In this notebook we'll be exploring the urban heat island effect by looking at the impact on surface temperature of roof color and street trees. We'll be replicating the process described here: http://urbanspatialanalysis.com/urban-heat-islands-street-trees-in-philadelphia/ but using Python tools rather than ESRI.

Extra packages: To run this notebook, you'll need the PyViz tools and a library of top of atmosphere calculations from rio-toa : pip install rio-toa

Data sources: This notebook uses Landsat data from Google Cloud Storage as well as some geographic data from OpenDataPhilly.

In [1]:
import intake
import xarray as xr
import pandas as pd
import numpy as np
import geopandas as gpd

import cartopy.crs as ccrs

import holoviews as hv
import hvplot.xarray
import hvplot.pandas

from geoviews.tile_sources import EsriImagery

hv.extension('bokeh')

Just some extra info about Landsat data:

In [2]:
band_info = pd.DataFrame([
        (1,  "Aerosol", " 0.43 - 0.45",    0.440,  "30",         "Coastal aerosol"),
        (2,  "Blue",    " 0.45 - 0.51",    0.480,  "30",         "Blue"),
        (3,  "Green",   " 0.53 - 0.59",    0.560,  "30",         "Green"),
        (4,  "Red",     " 0.64 - 0.67",    0.655,  "30",         "Red"),
        (5,  "NIR",     " 0.85 - 0.88",    0.865,  "30",         "Near Infrared (NIR)"),
        (6,  "SWIR1",   " 1.57 - 1.65",    1.610,  "30",         "Shortwave Infrared (SWIR) 1"),
        (7,  "SWIR2",   " 2.11 - 2.29",    2.200,  "30",         "Shortwave Infrared (SWIR) 2"),
        (8,  "Panc",    " 0.50 - 0.68",    0.590,  "15",         "Panchromatic"),
        (9,  "Cirrus",  " 1.36 - 1.38",    1.370,  "30",         "Cirrus"),
        (10, "TIRS1",   "10.60 - 11.19",   10.895, "100 * (30)", "Thermal Infrared (TIRS) 1"),
        (11, "TIRS2",   "11.50 - 12.51",   12.005, "100 * (30)", "Thermal Infrared (TIRS) 2")],
    columns=['Band', 'Name', 'Wavelength Range (µm)', 'Nominal Wavelength (µm)', 'Resolution (m)', 'Description']).set_index(["Band"])
band_info
Out[2]:
Name Wavelength Range (µm) Nominal Wavelength (µm) Resolution (m) Description
Band
1 Aerosol 0.43 - 0.45 0.440 30 Coastal aerosol
2 Blue 0.45 - 0.51 0.480 30 Blue
3 Green 0.53 - 0.59 0.560 30 Green
4 Red 0.64 - 0.67 0.655 30 Red
5 NIR 0.85 - 0.88 0.865 30 Near Infrared (NIR)
6 SWIR1 1.57 - 1.65 1.610 30 Shortwave Infrared (SWIR) 1
7 SWIR2 2.11 - 2.29 2.200 30 Shortwave Infrared (SWIR) 2
8 Panc 0.50 - 0.68 0.590 15 Panchromatic
9 Cirrus 1.36 - 1.38 1.370 30 Cirrus
10 TIRS1 10.60 - 11.19 10.895 100 * (30) Thermal Infrared (TIRS) 1
11 TIRS2 11.50 - 12.51 12.005 100 * (30) Thermal Infrared (TIRS) 2

Loading data

For this example, we will be using landsat data stored on Google Cloud Storage. Since these data are accessed via https, there is no guaranteed directory structure, so we will need to specify the url pointing to each file and then iterate over the files to create a concatenated dataset. We use jinja template notation in intake to pass parameters to the urlpath .

In [3]:
cat = intake.open_catalog('../catalog.yml')
list(cat)
Out[3]:
['l5',
 'l8',
 'google_landsat_band',
 'amazon_landsat_band',
 'fluxnet_daily',
 'fluxnet_metadata']

Let's take a look at what the google_landsat_band looks like:

yml
google_landsat_band:
    description: Landsat bands from Google Cloud Storage
    driver: rasterio
    parameters:
      path:
        description: landsat path
        type: int
      row:
        description: landsat row
        type: int
      product_id:
        description: landsat file id
        type: str
      band:
        description: band
        type: int
    args:
      urlpath: https://storage.googleapis.com/gcp-public-data-landsat/LC08/01/{{ '%03d' % path }}/{{ '%03d' % row }}/{{ product_id }}/{{ product_id }}_B{{ band }}.TIF
      chunks:
        band: 1
        x: 256
        y: 256

The following might feel a bit arbitrary, but we have chosen the path and row corresponding to the area over Philadelphia using the earth explorer . We have also found the id of the particular date of interest using the same tool. With these values in hand, we can access parts of each file on Google Cloud Storage.

In [4]:
path = 14
row = 32
product_id = 'LC08_L1TP_014032_20160727_20170222_01_T1'

For this workflow, we don't ever need to load in all of the data, just the metadata and some of the actual data. We do need to use "backoff" to avoid transient error conditions for https access for the metadata:

In [5]:
from random import random
from time import sleep

def get_band_with_exponential_backoff(path, row, product_id, band, 
                                      maximum_backoff=32):
    """
    Google Cloud Storage recommends using exponential backoff 
    when accessing the API. 
    
    https://cloud.google.com/storage/docs/exponential-backoff
    """
    n = backoff = 0
    while backoff < maximum_backoff:
        try:
            return cat.google_landsat_band(path=path, row=row, 
                                           product_id=product_id,
                                           band=band).to_dask()
        except:
            backoff = min(2**n + random(), maximum_backoff)
            sleep(backoff)
            n += 1
In [6]:
# Load all the bands
# bands = range(1, 12)

# OR skip band 8 for now since it is at a different resolution and transform
bands = [1, 2, 3, 4, 5, 6, 7, 9, 10, 11]

datasets = []
for band in bands:
    da = get_band_with_exponential_backoff(path=path, row=row, product_id=product_id, band=band)
    da = da.assign_coords(band=[band])
    datasets.append(da)
    
ds = xr.concat(datasets, 'band', compat='identical')
print(ds)
<xarray.DataArray (band: 10, y: 7871, x: 7741)>
dask.array<shape=(10, 7871, 7741), dtype=uint16, chunksize=(1, 256, 256)>
Coordinates:
  * y        (y) float64 4.582e+06 4.582e+06 4.582e+06 ... 4.346e+06 4.346e+06
  * x        (x) float64 3.954e+05 3.954e+05 3.955e+05 ... 6.276e+05 6.276e+05
  * band     (band) int64 1 2 3 4 5 6 7 9 10 11
Attributes:
    transform:   (30.0, 0.0, 395385.0, 0.0, -30.0, 4582215.0)
    crs:         +init=epsg:32618
    res:         (30.0, 30.0)
    is_tiled:    1
    nodatavals:  (nan,)

Loading in metadata regarding these particular Landsat images from the associated matlab.txt file.

In [7]:
def load_google_landsat_metadata(path, row, product_id):
    """Load Landsat metadata for path, row, product_id from Google Cloud Storage
    """
    def parse_type(x):
        try: 
            return eval(x)
        except:
            return x
    
    baseurl = 'https://storage.googleapis.com/gcp-public-data-landsat/LC08/01'
    suffix = f'{path:03d}/{row:03d}/{product_id}/{product_id}_MTL.txt'
    
    df = intake.open_csv(
        urlpath=f'{baseurl}/{suffix}',
        csv_kwargs={'sep': '=',
                    'header': None,
                    'names': ['index', 'value'],
                    'skiprows': 2,
                    'converters': {'index': (lambda x: x.strip()),
                                   'value': parse_type}}).read()
    metadata = df.set_index('index')['value']
    return metadata
In [8]:
metadata = load_google_landsat_metadata(path, row, product_id)
metadata.head()
Out[8]:
index
ORIGIN                Image courtesy of the U.S. Geological Survey
REQUEST_ID                                     0501702206266_00020
LANDSAT_SCENE_ID                             LC80140322016209LGN01
LANDSAT_PRODUCT_ID        LC08_L1TP_014032_20160727_20170222_01_T1
COLLECTION_NUMBER                                               01
Name: value, dtype: object

Sub-setting to area of interest

So far we haven't downloaded any band data. Since we know that we are interested in Philadelphia, we can just take a smaller square of data that covers the extents of the city. First we need to know the projection of the dataset:

In [9]:
ds.crs
Out[9]:
'+init=epsg:32618'

We'll convert that into something directly usable for later:

In [10]:
crs = ccrs.epsg(32618)

Now if we were just looking for one particular point we could use that point, converted to the coordinate system of the data, and then select the data nearest to it:

In [11]:
x_center, y_center = crs.transform_point(-75.1652, 39.9526, ccrs.PlateCarree())
nearest_to_center = ds.sel(x=x_center, y=y_center, method='nearest')
print(nearest_to_center.compute())

nearest_to_center.hvplot()
<xarray.DataArray (band: 10)>
array([13010, 12392, 11576, 11465, 12774, 11526,  9693,  5050, 30347, 26335],
      dtype=uint16)
Coordinates:
    y        float64 4.423e+06
    x        float64 4.859e+05
  * band     (band) int64 1 2 3 4 5 6 7 9 10 11
Attributes:
    transform:   (30.0, 0.0, 395385.0, 0.0, -30.0, 4582215.0)
    crs:         +init=epsg:32618
    res:         (30.0, 30.0)
    is_tiled:    1
    nodatavals:  (nan,)
Out[11]:

In this case, though, we are interested in a subset of data that covers that city of Philadelphia. So we need some geometry to specify the bounds of the city. We can get a GeoJSON of neighborhood data from OpenDataPhilly .

In [12]:
url = 'https://github.com/azavea/geo-data/raw/master/Neighborhoods_Philadelphia/Neighborhoods_Philadelphia.geojson'
geoms = gpd.read_file(url)

The approach shown below is somewhat simplistic in that we will use the map to iteratively select the area of interest by tweaking the bounds out from the central lat, lon of interest.

In [13]:
subset = ds.sel(x=slice(x_center-1.2e4, x_center+1.9e4), y=slice(y_center+2.1e4, y_center-1.2e4))
In [14]:
plot = subset.sel(band=3).hvplot(rasterize=True, crs=crs, height=500, width=600).redim(x="Longitude", y="Latitude")
plot.options(tools=[]) * geoms.hvplot(geo=True, alpha=.5, c='mapname')
Out[14]:

Now that we have looked at our map and ensured that the area covers the city, we will just take this chunk of the bands to use for further analyses:

In [15]:
ds = subset
ds
Out[15]:
<xarray.DataArray (band: 10, y: 1100, x: 1033)>
dask.array<shape=(10, 1100, 1033), dtype=uint16, chunksize=(1, 240, 199)>
Coordinates:
  * y        (y) float64 4.443e+06 4.443e+06 4.443e+06 ... 4.411e+06 4.411e+06
  * x        (x) float64 4.739e+05 4.739e+05 4.74e+05 ... 5.048e+05 5.049e+05
  * band     (band) int64 1 2 3 4 5 6 7 9 10 11
Attributes:
    transform:   (30.0, 0.0, 395385.0, 0.0, -30.0, 4582215.0)
    crs:         +init=epsg:32618
    res:         (30.0, 30.0)
    is_tiled:    1
    nodatavals:  (nan,)

Calculate NDVI

We'll calculate NDVI but we won't yet do any computations -- our bands are actually dask arrays, which allow for lazy computation.

In [16]:
NDVI = (ds.sel(band=5) - ds.sel(band=4)) / (ds.sel(band=5) + ds.sel(band=4))
NDVI = NDVI.where(NDVI < np.inf)
NDVI
Out[16]:
<xarray.DataArray (y: 1100, x: 1033)>
dask.array<shape=(1100, 1033), dtype=float64, chunksize=(240, 199)>
Coordinates:
  * y        (y) float64 4.443e+06 4.443e+06 4.443e+06 ... 4.411e+06 4.411e+06
  * x        (x) float64 4.739e+05 4.739e+05 4.74e+05 ... 5.048e+05 5.049e+05

In order to visualize NDVI, the data will need to be loaded and the NDVI computed. We can expect this to take some non-trivial amount of time (on the order of 20 sec on my machine).

In [17]:
p = NDVI.hvplot(datashade=True, x='x', y='y', crs=crs, height=500, width=500, cmap='viridis')

# We'll use a transparent rasterized version of the plot for hover text
p_hover = NDVI.hvplot(rasterize=True, x='x', y='y', crs=crs, height=500, width=500, alpha=0, colorbar=False)

p * p_hover
Out[17]:

Calculate land surface temperature

Given the NDVI calculated above, we can determine land surface temperature. For ease, we'll use some top of atmosphere calculations that have already been written up as Python functions as part of rasterio work in the rio_toa module. We'll also need to specify one more for transforming satellite temperature (brightness temp) to land surface temperature. For these calculations we'll use both Thermal Infrared bands - 10 and 11.

In [18]:
from rio_toa import brightness_temp, toa_utils
In [19]:
def land_surface_temp(BT, w, NDVI):
    """Calculate land surface temperature of Landsat 8
    
    temp = BT/1 + w * (BT /p) * ln(e)
    
    BT = At Satellite temperature (brightness)
    w = wavelength of emitted radiance (μm)
    
    where p = h * c / s (1.439e-2 mK)
    
    h = Planck's constant (Js)
    s = Boltzmann constant (J/K)
    c = velocity of light (m/s)
    """
    h = 6.626e-34
    s = 1.38e-23
    c = 2.998e8
    
    p = (h * c / s) * 1e6
    
    Pv = (NDVI - NDVI.min() / NDVI.max() - NDVI.min())**2
    e = 0.004 * Pv + 0.986
    
    return BT / 1 + w * (BT / p) * np.log(e)

Now we'll set up a helper function to retrieve all the parameters from the metadata and general Landsat info table, and calculate the land surface temperature for bands 10 and 11.

In [20]:
def land_surface_temp_for_band(band):
    # params from general Landsat info
    w = band_info.loc[band]['Nominal Wavelength (µm)']
    
    # params from specific Landsat data text file for these images
    ML = metadata[f'RADIANCE_MULT_BAND_{band}']
    AL = metadata[f'RADIANCE_ADD_BAND_{band}']
    K1 = metadata[f'K1_CONSTANT_BAND_{band}']
    K2 = metadata[f'K2_CONSTANT_BAND_{band}']
    
    at_satellite_temp = brightness_temp.brightness_temp(ds.sel(band=band).values, ML, AL, K1, K2)
    
    return land_surface_temp(at_satellite_temp, w, NDVI)
In [21]:
band = 10
temp_10 = land_surface_temp_for_band(band).compute()
temp_10_f = toa_utils.temp_rescale(temp_10, 'F')
band_10=ds.sel(band=band).copy(data=temp_10_f).hvplot(rasterize=True, x='x', y='y', cmap='fire_r', 
                                                      crs=crs, height=500, width=450)

band = 11
temp_11 = land_surface_temp_for_band(band).compute()
temp_11_f = toa_utils.temp_rescale(temp_11, 'F')
band_11=ds.sel(band=band).copy(data=temp_11_f).hvplot(rasterize=True, x='x', y='y', cmap='fire_r', 
                                                      crs=crs, height=500, width=450)

Compare the results from the two different bands, noticing that the colorbars are different.

In [22]:
band_10.relabel('Band 10') + band_11.relabel('Band 11')
Out[22]:

We'll take the mean of the calculated land surface temperature for each of the bands and mimic the colormap used in the project that we are duplicating.

In [23]:
mean_temp = (temp_10 + temp_11) / 2
mean_temp_f = toa_utils.temp_rescale(mean_temp, 'F')

data = ds[0].copy(data=mean_temp_f)
p = (data
    .hvplot(rasterize=True, x='x', y='y', crs=crs, height=500, width=600, cmap='rainbow', alpha=.5, legend=False)
    .relabel('Mean Surface Temp (F)'))

p * EsriImagery 
Out[23]:

Notice how the hot spots are located over warehouse roofs and parking lots. This becomes even more visible when just the temperatures greater than 90F are displayed. To show this, we'll make a special colormap that just includes high intensity reds that are found at the top of the fire_r colormap.

In [24]:
import colorcet as cc

special_cmap = cc.fire[::-1][90:]
In [25]:
data = ds[0].copy(data=mean_temp_f.where(mean_temp_f > 90))
thresholded_temp_p = (data
    .hvplot(x='x', y='y', cmap=special_cmap, crs=crs, height=500, width=450, 
            colorbar=False, legend=False, rasterize=True)
    .relabel('Mean Temp (F) > 90')
    .redim(value='Temperature (F)'))

thresholded_temp_p + thresholded_temp_p.options(alpha=.3) * EsriImagery
Out[25]:

Adding in the Street Tree data

OpenDataPhilly released an inventory of all the street trees in the city. Street trees are trees that are planted along streets, not those in parks and private property. The original analysis considered these 100,000 points too many to plot, but that's nothing to datashader , which is happy with billions of points even on a laptop.

It is hypothesized that where there are more trees the land surface temperature will be less extreme. To explore this, we will overlay street trees with the thresholded land surface temperature:

In [26]:
url = 'http://data.phl.opendata.arcgis.com/datasets/957f032f9c874327a1ad800abd887d17_0.geojson'
trees = gpd.read_file(url)
import geoviews as gv
gv.Points.datatype = ['geodataframe', 'dictionary'] # Temporary geoviews workaround
trees.tail()
Out[26]:
OBJECTID SPECIES STATUS DBH geometry
111988 111989 None None None POINT (-75.15905471478176 39.93696759022805)
111989 111990 None None None POINT (-75.1596428445722 39.93715882988963)
111990 111991 None None None POINT (-75.15979123336561 39.93719121929897)
111991 111992 None None None POINT (-75.14504238099485 39.94169632627076)
111992 111993 None None None POINT (-75.14451560227096 39.94174213633205)
In [27]:
tree_p = trees.hvplot(geo=True, datashade=True, height=500, width=400, legend=False, dynspread=30).relabel('Street Tree Density')

thresholded_temp_p.options(alpha=.5) * tree_p.options(alpha=.5)
Out[27]: